This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
Execute chunks of code by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.
Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.
When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).
library(dplyr)
Attaching package: ‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
library(tidyr)
library(ggplot2)
library(rgdal)
Loading required package: sp
rgdal: version: 1.4-6, (SVN revision 841)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/gdal
GDAL binary built with GEOS: FALSE
Loaded PROJ.4 runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/proj
Linking to sp version: 1.3-1
library(sp)
library(raster)
Attaching package: ‘raster’
The following object is masked from ‘package:tidyr’:
extract
The following object is masked from ‘package:dplyr’:
select
library(sf)
Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
list.files()
[1] "data" "r_geospatial_notebook.nb.html"
[3] "r_geospatial_notebook.Rmd" "README.md"
GDALinfo("./data/HARV_dsmCrop.tif")
rows 1367
columns 1697
bands 1
lower left origin.x 731453
lower left origin.y 4712471
res.x 1
res.y 1
ysign -1
oblique.x 0
oblique.y 0
driver GTiff
projection +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
file ./data/HARV_dsmCrop.tif
apparent band summary:
GDType hasNoDataValue NoDataValue blockSize1 blockSize2
1 Float64 TRUE -9999 1 1697
apparent band statistics:
Bmin Bmax Bmean Bsd
1 305.07 416.07 359.8531 17.83169
Metadata:
AREA_OR_POINT=Area
dsm <- raster("./data/HARV_dsmCrop.tif")
dsm
class : RasterLayer
dimensions : 1367, 1697, 2319799 (nrow, ncol, ncell)
resolution : 1, 1 (x, y)
extent : 731453, 733150, 4712471, 4713838 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : /Users/francesdavenport/Desktop/r-geospatial-100720-master/data/HARV_dsmCrop.tif
names : HARV_dsmCrop
values : 305.07, 416.07 (min, max)
class(dsm)
[1] "RasterLayer"
attr(,"package")
[1] "raster"
plot(dsm)
summary(dsm)
summary is an estimate based on a sample of 1e+05 cells (4.31% of all cells)
HARV_dsmCrop
Min. 305.56
1st Qu. 345.67
Median 359.70
3rd Qu. 374.24
Max. 415.40
NA's 0.00
cellStats(dsm, stat='mean')
[1] 359.8531
dsm[1:8, 1:8]
[1] 408.77 408.23 406.53 406.61 409.06 409.45 413.97 415.77 407.05 406.62 404.98
[12] 405.01 407.44 407.17 413.21 414.12 407.06 406.03 403.55 404.30 408.55 411.42
[23] 409.59 409.82 405.70 406.32 404.16 403.81 404.19 403.71 405.83 408.20 404.57
[34] 412.16 396.06 392.26 401.96 402.64 397.24 396.84 407.40 413.08 412.37 407.52
[45] 407.39 405.97 399.11 403.42 411.96 413.86 413.08 414.37 410.71 406.57 399.31
[56] 402.24 412.48 414.86 414.68 415.77 414.16 414.14 401.23 413.67
crs(dsm)
CRS arguments:
+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
+towgs84=0,0,0
dsm_df <- as.data.frame(dsm, xy=TRUE)
head(dsm_df)
x y HARV_dsmCrop
1 731453.5 4713838 408.77
2 731454.5 4713838 408.23
3 731455.5 4713838 406.53
4 731456.5 4713838 406.61
5 731457.5 4713838 409.06
6 731458.5 4713838 409.45
nrow(dsm)
[1] 1367
ncol(dsm)
[1] 1697
nrow(dsm_df)
[1] 2319799
1367*1697
[1] 2319799
ggplot() +
geom_raster(data = dsm_df, aes(x=x, y=y, fill=HARV_dsmCrop)) +
scale_fill_viridis_c()
ggplot() +
geom_histogram(data = dsm_df, aes(HARV_dsmCrop), bins=40)
dsm_df <- dsm_df %>%
mutate(elevation_bin = cut(HARV_dsmCrop, breaks = c(300, 340, 380, 420)))
ggplot() +
geom_bar(data = dsm_df, aes(elevation_bin))
dtm <- raster("./data/HARV_dtmCrop.tif")
plot(dtm)
canopy_height <- dsm - dtm
plot(canopy_height)
mean_dtm <- cellStats(dtm, stat='mean')
mean_dtm
[1] 344.8979
dtm_anom <- dtm - mean_dtm
plot(dtm_anom)
writeRaster(canopy_height, "./data/canopy_height.tif", overwrite=TRUE)
harv_rgb <- stack("./data/HARV_RGB_Ortho.tif")
harv_rgb
class : RasterStack
dimensions : 2317, 3073, 7120141, 3 (nrow, ncol, ncell, nlayers)
resolution : 0.25, 0.25 (x, y)
extent : 731998.5, 732766.8, 4712956, 4713536 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
names : HARV_RGB_Ortho.1, HARV_RGB_Ortho.2, HARV_RGB_Ortho.3
min values : 0, 0, 0
max values : 255, 255, 255
plot(harv_rgb)
plotRGB(harv_rgb, r = 1, g = 2, b = 3)
aoi_sp <- shapefile("./data/HarClip_UTMZ18.shp")
aoi_sp
class : SpatialPolygonsDataFrame
features : 1
extent : 732128, 732251.1, 4713209, 4713359 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
variables : 1
names : id
value : 1
plot(aoi_sp)
aoi_sp[1,]
class : SpatialPolygonsDataFrame
features : 1
extent : 732128, 732251.1, 4713209, 4713359 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
variables : 1
names : id
value : 1
aoi <- st_read("./data/HarClip_UTMZ18.shp")
Reading layer `HarClip_UTMZ18' from data source `/Users/francesdavenport/Desktop/r-geospatial-100720-master/data/HarClip_UTMZ18.shp' using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 1 field
geometry type: POLYGON
dimension: XY
bbox: xmin: 732128 ymin: 4713209 xmax: 732251.1 ymax: 4713359
epsg (SRID): 32618
proj4string: +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
plot(aoi)
head(aoi)
Simple feature collection with 1 feature and 1 field
geometry type: POLYGON
dimension: XY
bbox: xmin: 732128 ymin: 4713209 xmax: 732251.1 ymax: 4713359
epsg (SRID): 32618
proj4string: +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
id geometry
1 1 POLYGON ((732128 4713359, 7...
ggplot() +
geom_sf(data = aoi)
ne_state <- st_read("./data/Boundary-US-State-NEast.shp")
Reading layer `Boundary-US-State-NEast' from data source `/Users/francesdavenport/Desktop/r-geospatial-100720-master/data/Boundary-US-State-NEast.shp' using driver `ESRI Shapefile'
Simple feature collection with 12 features and 9 fields
geometry type: MULTIPOLYGON
dimension: XYZ
bbox: xmin: -80.51989 ymin: 37.91685 xmax: -66.9499 ymax: 47.45716
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
plot(ne_state)
names(ne_state)
[1] "STATEFP" "STATENS" "AFFGEOID" "GEOID" "STUSPS" "NAME" "LSAD"
[8] "ALAND" "AWATER" "geometry"
head(ne_state)
Simple feature collection with 6 features and 9 fields
geometry type: MULTIPOLYGON
dimension: XYZ
bbox: xmin: -79.76212 ymin: 37.91685 xmax: -66.9499 ymax: 47.45716
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
STATEFP STATENS AFFGEOID GEOID STUSPS NAME LSAD ALAND
1 11 01702382 0400000US11 11 DC District of Columbia 00 158350578
2 24 01714934 0400000US24 24 MD Maryland 00 25147575220
3 36 01779796 0400000US36 36 NY New York 00 122054577774
4 09 01779780 0400000US09 09 CT Connecticut 00 12542396439
5 23 01779787 0400000US23 23 ME Maine 00 79885244456
6 25 00606926 0400000US25 25 MA Massachusetts 00 20203936287
AWATER geometry
1 18633500 MULTIPOLYGON Z (((-77.11976...
2 6983455225 MULTIPOLYGON Z (((-76.04621...
3 19242052501 MULTIPOLYGON Z (((-72.01893...
4 1814978794 MULTIPOLYGON Z (((-73.69594...
5 11749091526 MULTIPOLYGON Z (((-68.92401...
6 7131805598 MULTIPOLYGON Z (((-70.27553...
ggplot() +
geom_sf(data = ne_state, aes(fill = NAME))
ggplot() +
geom_sf(data = ne_state, aes(fill = AWATER))
mass <- filter(ne_state, NAME == "Massachusetts")
ggplot() +
geom_sf(data = mass)
mass
Simple feature collection with 1 feature and 9 fields
geometry type: MULTIPOLYGON
dimension: XYZ
bbox: xmin: -73.50814 ymin: 41.23796 xmax: -69.92826 ymax: 42.88459
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
STATEFP STATENS AFFGEOID GEOID STUSPS NAME LSAD ALAND AWATER
1 25 00606926 0400000US25 25 MA Massachusetts 00 20203936287 7131805598
geometry
1 MULTIPOLYGON Z (((-70.27553...
roads <- st_read("./data/HARV_roads.shp")
Reading layer `HARV_roads' from data source `/Users/francesdavenport/Desktop/r-geospatial-100720-master/data/HARV_roads.shp' using driver `ESRI Shapefile'
Simple feature collection with 13 features and 15 fields
geometry type: MULTILINESTRING
dimension: XY
bbox: xmin: 730741.2 ymin: 4711942 xmax: 733295.5 ymax: 4714260
epsg (SRID): 32618
proj4string: +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
names(roads)
[1] "OBJECTID_1" "OBJECTID" "TYPE" "NOTES" "MISCNOTES" "RULEID"
[7] "MAPLABEL" "SHAPE_LENG" "LABEL" "BIKEHORSE" "RESVEHICLE" "RECMAP"
[13] "Shape_Le_1" "ResVehic_1" "BicyclesHo" "geometry"
head(roads)
Simple feature collection with 6 features and 15 fields
geometry type: MULTILINESTRING
dimension: XY
bbox: xmin: 730741.2 ymin: 4712685 xmax: 732232.3 ymax: 4713726
epsg (SRID): 32618
proj4string: +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
OBJECTID_1 OBJECTID TYPE NOTES MISCNOTES RULEID MAPLABEL
1 14 48 woods road Locust Opening Rd <NA> 5 Locust Opening Rd
2 40 91 footpath <NA> <NA> 6 <NA>
3 41 106 footpath <NA> <NA> 6 <NA>
4 211 279 stone wall <NA> <NA> 1 <NA>
5 212 280 stone wall <NA> <NA> 1 <NA>
6 213 281 stone wall <NA> <NA> 1 <NA>
SHAPE_LENG LABEL BIKEHORSE RESVEHICLE RECMAP Shape_Le_1
1 1297.35706 Locust Opening Rd Y R1 Y 1297.10617
2 146.29984 <NA> Y R1 Y 146.29983
3 676.71804 <NA> Y R2 Y 676.71807
4 231.78957 <NA> <NA> <NA> <NA> 231.78962
5 45.50864 <NA> <NA> <NA> <NA> 45.50859
6 198.39043 <NA> <NA> <NA> <NA> 198.39041
ResVehic_1 BicyclesHo
1 R1 - All Research Vehicles Allowed Bicycles and Horses Allowed
2 R1 - All Research Vehicles Allowed Bicycles and Horses Allowed
3 R2 - 4WD/High Clearance Vehicles Only Bicycles and Horses Allowed
4 <NA> <NA>
5 <NA> <NA>
6 <NA> <NA>
geometry
1 MULTILINESTRING ((730819.2 ...
2 MULTILINESTRING ((732040.2 ...
3 MULTILINESTRING ((732057 47...
4 MULTILINESTRING ((731903.6 ...
5 MULTILINESTRING ((732039.1 ...
6 MULTILINESTRING ((732056.2 ...
roads[1,]
Simple feature collection with 1 feature and 15 fields
geometry type: MULTILINESTRING
dimension: XY
bbox: xmin: 730741.2 ymin: 4712685 xmax: 731960.4 ymax: 4713205
epsg (SRID): 32618
proj4string: +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
OBJECTID_1 OBJECTID TYPE NOTES MISCNOTES RULEID MAPLABEL
1 14 48 woods road Locust Opening Rd <NA> 5 Locust Opening Rd
SHAPE_LENG LABEL BIKEHORSE RESVEHICLE RECMAP Shape_Le_1
1 1297.357 Locust Opening Rd Y R1 Y 1297.106
ResVehic_1 BicyclesHo
1 R1 - All Research Vehicles Allowed Bicycles and Horses Allowed
geometry
1 MULTILINESTRING ((730819.2 ...
roads$TYPE
[1] woods road footpath footpath stone wall stone wall stone wall stone wall
[8] stone wall stone wall boardwalk woods road woods road woods road
Levels: boardwalk footpath stone wall woods road
roads[,"LABEL"]
Simple feature collection with 13 features and 1 field
geometry type: MULTILINESTRING
dimension: XY
bbox: xmin: 730741.2 ymin: 4711942 xmax: 733295.5 ymax: 4714260
epsg (SRID): 32618
proj4string: +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
First 10 features:
LABEL geometry
1 Locust Opening Rd MULTILINESTRING ((730819.2 ...
2 <NA> MULTILINESTRING ((732040.2 ...
3 <NA> MULTILINESTRING ((732057 47...
4 <NA> MULTILINESTRING ((731903.6 ...
5 <NA> MULTILINESTRING ((732039.1 ...
6 <NA> MULTILINESTRING ((732056.2 ...
7 <NA> MULTILINESTRING ((731964 47...
8 <NA> MULTILINESTRING ((732105.2 ...
9 <NA> MULTILINESTRING ((732222.9 ...
10 <NA> MULTILINESTRING ((732153.8 ...
st_crs(roads)
Coordinate Reference System:
EPSG: 32618
proj4string: "+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs"
st_geometry(roads)
Geometry set for 13 features
geometry type: MULTILINESTRING
dimension: XY
bbox: xmin: 730741.2 ymin: 4711942 xmax: 733295.5 ymax: 4714260
epsg (SRID): 32618
proj4string: +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
First 5 geometries:
MULTILINESTRING ((730819.2 4712714, 730812.1 47...
MULTILINESTRING ((732040.2 4713241, 732033.7 47...
MULTILINESTRING ((732057 4713726, 732053.8 4713...
MULTILINESTRING ((731903.6 4713197, 731882.6 47...
MULTILINESTRING ((732039.1 4713299, 732051.9 47...
st_coordinates(roads)
X Y L1 L2
[1,] 730819.2 4712714 1 1
[2,] 730812.1 4712714 1 1
[3,] 730809.8 4712714 1 1
[4,] 730808.7 4712713 1 1
[5,] 730807.1 4712712 1 1
[6,] 730805.4 4712710 1 1
[7,] 730798.7 4712700 1 1
[8,] 730797.2 4712699 1 1
[9,] 730795.7 4712697 1 1
[10,] 730793.1 4712695 1 1
[11,] 730790.4 4712693 1 1
[12,] 730787.6 4712692 1 1
[13,] 730784.8 4712690 1 1
[14,] 730781.7 4712689 1 1
[15,] 730778.6 4712688 1 1
[16,] 730775.4 4712687 1 1
[17,] 730772.0 4712686 1 1
[18,] 730768.6 4712686 1 1
[19,] 730765.0 4712685 1 1
[20,] 730757.5 4712685 1 1
[21,] 730749.6 4712686 1 1
[22,] 730741.2 4712687 1 1
[23,] 731034.8 4712792 2 1
[24,] 731042.7 4712794 2 1
[25,] 731050.4 4712795 2 1
[26,] 731058.0 4712798 2 1
[27,] 731065.4 4712800 2 1
[28,] 731075.9 4712804 2 1
[29,] 731083.0 4712807 2 1
[30,] 731089.2 4712810 2 1
[31,] 731094.4 4712814 2 1
[32,] 731097.9 4712817 2 1
[33,] 731101.6 4712821 2 1
[34,] 731105.4 4712825 2 1
[35,] 731111.4 4712833 2 1
[36,] 731115.7 4712838 2 1
[37,] 731120.2 4712845 2 1
[38,] 731133.2 4712865 2 1
[39,] 731135.8 4712870 2 1
[40,] 731143.3 4712883 2 1
[41,] 731146.0 4712888 2 1
[42,] 731148.6 4712891 2 1
[43,] 731152.1 4712895 2 1
[44,] 731154.9 4712898 2 1
[45,] 731157.0 4712900 2 1
[46,] 731159.4 4712902 2 1
[47,] 731162.2 4712903 2 1
[48,] 731165.2 4712904 2 1
[49,] 731171.6 4712906 2 1
[50,] 731184.8 4712909 2 1
[51,] 731190.2 4712910 2 1
[52,] 731195.5 4712912 2 1
[53,] 731199.9 4712914 2 1
[54,] 731204.3 4712917 2 1
[55,] 731208.9 4712920 2 1
[56,] 731212.8 4712924 2 1
[57,] 731216.7 4712928 2 1
[58,] 731220.8 4712932 2 1
[59,] 731232.3 4712946 2 1
[60,] 731236.2 4712950 2 1
[61,] 731240.0 4712954 2 1
[62,] 731247.2 4712960 2 1
[63,] 731249.7 4712962 2 1
[64,] 731252.0 4712963 2 1
[65,] 731255.0 4712964 2 1
[66,] 731258.3 4712966 2 1
[67,] 731261.9 4712967 2 1
[68,] 731265.8 4712968 2 1
[69,] 731269.1 4712968 2 1
[70,] 731272.9 4712968 2 1
[71,] 731285.3 4712967 2 1
[72,] 731289.8 4712967 2 1
[73,] 731298.3 4712968 2 1
[74,] 731305.6 4712968 2 1
[75,] 731312.2 4712969 2 1
[76,] 731318.3 4712970 2 1
[77,] 731325.9 4712971 2 1
[78,] 731332.4 4712973 2 1
[79,] 731338.2 4712974 2 1
[80,] 731343.5 4712976 2 1
[81,] 731350.0 4712979 2 1
[82,] 731355.6 4712981 2 1
[83,] 731360.6 4712984 2 1
[84,] 731365.0 4712987 2 1
[85,] 731367.5 4712989 2 1
[86,] 731370.3 4712991 2 1
[87,] 731378.8 4713000 2 1
[88,] 731382.8 4713004 2 1
[89,] 731393.1 4713016 2 1
[90,] 731396.5 4713019 2 1
[91,] 731399.8 4713022 2 1
[92,] 731406.3 4713026 2 1
[93,] 731414.7 4713030 2 1
[94,] 731417.9 4713032 2 1
[95,] 731421.6 4713033 2 1
[96,] 731434.2 4713036 2 1
[97,] 731438.6 4713038 2 1
[98,] 731442.7 4713039 2 1
[99,] 731446.1 4713041 2 1
[100,] 731453.4 4713046 2 1
[101,] 731467.1 4713056 2 1
[102,] 731476.4 4713062 2 1
[103,] 731485.1 4713067 2 1
[104,] 731493.3 4713070 2 1
[105,] 731496.8 4713072 2 1
[106,] 731500.8 4713073 2 1
[107,] 731514.3 4713076 2 1
[108,] 731519.1 4713077 2 1
[109,] 731523.6 4713078 2 1
[110,] 731527.5 4713080 2 1
[111,] 731532.7 4713083 2 1
[112,] 731538.2 4713087 2 1
[113,] 731547.9 4713095 2 1
[114,] 731552.9 4713098 2 1
[115,] 731559.9 4713102 2 1
[116,] 731567.8 4713106 2 1
[117,] 731577.0 4713110 2 1
[118,] 731586.9 4713113 2 1
[119,] 731604.9 4713118 2 1
[120,] 731614.1 4713121 2 1
[121,] 731619.3 4713122 2 1
[122,] 731623.3 4713123 2 1
[123,] 731633.7 4713123 2 1
[124,] 731644.4 4713124 2 1
[125,] 731652.4 4713125 2 1
[126,] 731659.5 4713127 2 1
[127,] 731665.8 4713129 2 1
[128,] 731674.2 4713132 2 1
[129,] 731677.4 4713133 2 1
[130,] 731680.2 4713135 2 1
[131,] 731682.4 4713137 2 1
[132,] 731687.1 4713141 2 1
[133,] 731689.5 4713142 2 1
[134,] 731692.5 4713144 2 1
[135,] 731696.4 4713145 2 1
[136,] 731706.3 4713148 2 1
[137,] 731711.4 4713150 2 1
[138,] 731719.4 4713154 2 1
[139,] 731725.6 4713157 2 1
[140,] 731730.1 4713159 2 1
[141,] 731734.6 4713162 2 1
[142,] 731746.1 4713171 2 1
[143,] 731758.1 4713178 2 1
[144,] 731768.9 4713185 2 1
[145,] 731779.2 4713190 2 1
[146,] 731789.1 4713195 2 1
[147,] 731796.6 4713197 2 1
[148,] 731803.6 4713200 2 1
[149,] 731810.9 4713202 2 1
[150,] 731815.8 4713203 2 1
[151,] 731820.9 4713204 2 1
[152,] 731827.8 4713205 2 1
[153,] 731834.5 4713205 2 1
[154,] 731841.9 4713204 2 1
[155,] 731848.7 4713204 2 1
[156,] 731852.9 4713204 2 1
[157,] 731857.3 4713202 2 1
[158,] 731865.1 4713199 2 1
[159,] 731880.8 4713194 2 1
[160,] 731891.1 4713190 2 1
[161,] 731900.3 4713186 2 1
[162,] 731908.3 4713182 2 1
[163,] 731915.0 4713178 2 1
[164,] 731920.8 4713174 2 1
[165,] 731926.0 4713170 2 1
[166,] 731930.7 4713166 2 1
[167,] 731936.2 4713162 2 1
[168,] 731939.8 4713158 2 1
[169,] 731942.9 4713155 2 1
[170,] 731945.3 4713152 2 1
[171,] 731947.9 4713148 2 1
[172,] 731950.2 4713143 2 1
[173,] 731952.3 4713138 2 1
[174,] 731954.5 4713131 2 1
[175,] 731955.7 4713127 2 1
[176,] 731956.7 4713122 2 1
[177,] 731958.7 4713105 2 1
[178,] 731959.9 4713094 2 1
[179,] 731960.4 4713086 2 1
[180,] 731960.3 4713079 2 1
[181,] 731959.7 4713069 2 1
[182,] 731959.2 4713062 2 1
[183,] 731957.4 4713047 2 1
[184,] 731956.9 4713040 2 1
[185,] 731956.5 4713030 2 1
[186,] 731956.4 4713025 2 1
[187,] 731956.8 4713019 2 1
[188,] 731957.5 4713015 2 1
[189,] 731958.6 4713010 2 1
[190,] 731960.1 4713006 2 1
[191,] 732040.2 4713241 1 2
[192,] 732033.7 4713235 1 2
[193,] 732022.1 4713230 1 2
[194,] 732013.8 4713222 1 2
[195,] 732009.0 4713212 1 2
[196,] 732004.9 4713200 1 2
[197,] 732002.9 4713187 1 2
[198,] 731999.3 4713175 1 2
[199,] 731993.5 4713164 1 2
[200,] 731983.4 4713149 1 2
[201,] 731971.3 4713141 1 2
[202,] 731954.5 4713131 1 2
[203,] 732057.0 4713726 1 3
[204,] 732053.8 4713719 1 3
[205,] 732053.9 4713711 1 3
[206,] 732054.9 4713703 1 3
[207,] 732058.7 4713694 1 3
[208,] 732064.4 4713680 1 3
[209,] 732071.0 4713673 1 3
[210,] 732076.7 4713659 1 3
[211,] 732078.8 4713644 1 3
[212,] 732079.9 4713635 1 3
[213,] 732081.8 4713626 1 3
[214,] 732102.2 4713617 1 3
[215,] 732143.9 4713602 1 3
[216,] 732167.2 4713590 1 3
[217,] 732200.2 4713560 1 3
[218,] 732221.3 4713531 1 3
[219,] 732225.2 4713499 1 3
[220,] 732219.1 4713477 1 3
[221,] 732213.8 4713460 1 3
[222,] 732211.2 4713446 1 3
[223,] 732226.9 4713416 1 3
[224,] 732232.3 4713401 1 3
[225,] 732229.3 4713387 1 3
[226,] 732206.2 4713354 1 3
[227,] 732157.0 4713309 1 3
[228,] 732138.3 4713290 1 3
[229,] 732124.1 4713282 1 3
[230,] 732108.8 4713278 1 3
[231,] 732095.5 4713275 1 3
[232,] 732079.2 4713259 1 3
[233,] 732061.1 4713245 1 3
[234,] 732046.4 4713242 1 3
[235,] 732040.2 4713241 1 3
[236,] 731903.6 4713197 1 4
[237,] 731882.6 4713287 1 4
[238,] 732021.8 4713296 1 4
[239,] 732039.1 4713299 1 5
[240,] 732051.9 4713255 1 5
[241,] 732056.2 4713235 1 6
[242,] 732105.2 4713043 1 6
[243,] 731964.0 4713019 1 7
[244,] 732105.2 4713043 1 7
[245,] 732105.2 4713043 1 8
[246,] 732194.2 4713059 1 8
[247,] 732222.9 4713063 1 9
[248,] 732258.2 4713070 1 9
[249,] 732153.8 4713305 1 10
[250,] 732159.8 4713303 1 10
[ reached getOption("max.print") -- omitted 109 rows ]
tower <- st_read("./data/HARVtower_UTM18N.shp")
Reading layer `HARVtower_UTM18N' from data source `/Users/francesdavenport/Desktop/r-geospatial-100720-master/data/HARVtower_UTM18N.shp' using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 14 fields
geometry type: POINT
dimension: XY
bbox: xmin: 732183.2 ymin: 4713265 xmax: 732183.2 ymax: 4713265
epsg (SRID): 32618
proj4string: +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
ggplot() +
geom_sf(data = aoi, fill = "green", alpha = 0.5) +
geom_sf(data = roads) +
geom_sf(data = tower, col = "red", size = 10, shape = "*")
ggplot() +
geom_sf(data = roads, aes(col = LABEL, linetype = TYPE))
tower_buffer <- st_buffer(tower, 200)
ggplot()+
geom_sf(data = tower) +
geom_sf(data = tower_buffer, fill = "blue", alpha = 0.3)
crs(canopy_height)
CRS arguments:
+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
st_crs(tower_buffer)
Coordinate Reference System:
EPSG: 32618
proj4string: "+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs"
canopy_masked <- mask(canopy_height, tower_buffer, inverse=TRUE)
plot(canopy_masked)
st_bbox(tower_buffer)
xmin ymin xmax ymax
731983.2 4713065.0 732383.2 4713465.0
canopy_masked <- crop(canopy_masked, extent(731900, 732500, 4713000, 4713500))
plot(canopy_masked)
cellStats(canopy_masked, stat = "mean")
[1] 19.36194
canopy_extracted_values <- raster::extract(canopy_height, tower_buffer)
canopy_extracted_values <- unlist(canopy_extracted_values)
hist(canopy_extracted_values)
ndvi_files <- list.files("./data/NDVI", full.names=TRUE, pattern = ".tif$")
ndvi_files
[1] "./data/NDVI/005_HARV_ndvi_crop.tif" "./data/NDVI/037_HARV_ndvi_crop.tif"
[3] "./data/NDVI/085_HARV_ndvi_crop.tif" "./data/NDVI/133_HARV_ndvi_crop.tif"
[5] "./data/NDVI/181_HARV_ndvi_crop.tif" "./data/NDVI/197_HARV_ndvi_crop.tif"
[7] "./data/NDVI/213_HARV_ndvi_crop.tif" "./data/NDVI/229_HARV_ndvi_crop.tif"
[9] "./data/NDVI/245_HARV_ndvi_crop.tif" "./data/NDVI/261_HARV_ndvi_crop.tif"
[11] "./data/NDVI/277_HARV_ndvi_crop.tif" "./data/NDVI/293_HARV_ndvi_crop.tif"
[13] "./data/NDVI/309_HARV_ndvi_crop.tif"
ndvi <- stack(ndvi_files)
ndvi
class : RasterStack
dimensions : 5, 4, 20, 13 (nrow, ncol, ncell, nlayers)
resolution : 30, 30 (x, y)
extent : 239415, 239535, 4714215, 4714365 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=19 +ellps=WGS84 +units=m +no_defs
names : X005_HARV_ndvi_crop, X037_HARV_ndvi_crop, X085_HARV_ndvi_crop, X133_HARV_ndvi_crop, X181_HARV_ndvi_crop, X197_HARV_ndvi_crop, X213_HARV_ndvi_crop, X229_HARV_ndvi_crop, X245_HARV_ndvi_crop, X261_HARV_ndvi_crop, X277_HARV_ndvi_crop, X293_HARV_ndvi_crop, X309_HARV_ndvi_crop
min values : 2732, 1534, 1917, 5669, 8685, 8768, 8633, 8686, 8297, 7491, 277, 484, 4396
max values : 5545, 3736, 3600, 6394, 8903, 9140, 8945, 8967, 8651, 8607, 366, 661, 6359
ndvi_Z18 <- projectRaster(ndvi, crs = crs(dtm))
no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf
ndvi_Z18
class : RasterBrick
dimensions : 5, 4, 20, 13 (nrow, ncol, ncell, nlayers)
resolution : 29.9, 29.9 (x, y)
extent : 732135.1, 732254.7, 4713216, 4713365 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : X005_HARV_ndvi_crop, X037_HARV_ndvi_crop, X085_HARV_ndvi_crop, X133_HARV_ndvi_crop, X181_HARV_ndvi_crop, X197_HARV_ndvi_crop, X213_HARV_ndvi_crop, X229_HARV_ndvi_crop, X245_HARV_ndvi_crop, X261_HARV_ndvi_crop, X277_HARV_ndvi_crop, X293_HARV_ndvi_crop, X309_HARV_ndvi_crop
min values : 2746.0677, 1649.4942, 2037.4145, 5752.6465, 8673.7318, 8813.5836, 8583.8179, 8711.7698, 8278.1903, 7542.8330, 295.9601, 461.3890, 4732.6922
max values : 5715.4455, 3872.7388, 3724.4348, 6445.5631, 8869.5036, 9138.6236, 8963.2645, 8944.5587, 8610.9731, 8405.7861, 368.3793, 630.1992, 6470.4486
plot(ndvi_Z18, 1)
cellStats(ndvi_Z18, stat="mean")
X005_HARV_ndvi_crop X037_HARV_ndvi_crop X085_HARV_ndvi_crop X133_HARV_ndvi_crop
3709.2514 2507.5816 2571.5682 6004.5853
X181_HARV_ndvi_crop X197_HARV_ndvi_crop X213_HARV_ndvi_crop X229_HARV_ndvi_crop
8780.7595 8941.2789 8783.9587 8818.7488
X245_HARV_ndvi_crop X261_HARV_ndvi_crop X277_HARV_ndvi_crop X293_HARV_ndvi_crop
8503.5736 7955.7723 332.0345 561.2631
X309_HARV_ndvi_crop
5457.6267
tower_ndvi <- raster::extract(ndvi_Z18, tower)
tower_ndvi
X005_HARV_ndvi_crop X037_HARV_ndvi_crop X085_HARV_ndvi_crop X133_HARV_ndvi_crop
[1,] 3253.704 1888.675 2631.495 5887.102
X181_HARV_ndvi_crop X197_HARV_ndvi_crop X213_HARV_ndvi_crop X229_HARV_ndvi_crop
[1,] 8793.096 8942.353 8867.162 8856.346
X245_HARV_ndvi_crop X261_HARV_ndvi_crop X277_HARV_ndvi_crop X293_HARV_ndvi_crop
[1,] 8583.951 7804.457 331.571 563.1186
X309_HARV_ndvi_crop
[1,] 5366.162
ndvi_df <- as.data.frame(tower_ndvi[1,])
ndvi_df
tower_ndvi[1, ]
X005_HARV_ndvi_crop 3253.7035
X037_HARV_ndvi_crop 1888.6754
X085_HARV_ndvi_crop 2631.4945
X133_HARV_ndvi_crop 5887.1022
X181_HARV_ndvi_crop 8793.0962
X197_HARV_ndvi_crop 8942.3530
X213_HARV_ndvi_crop 8867.1619
X229_HARV_ndvi_crop 8856.3461
X245_HARV_ndvi_crop 8583.9509
X261_HARV_ndvi_crop 7804.4574
X277_HARV_ndvi_crop 331.5710
X293_HARV_ndvi_crop 563.1186
X309_HARV_ndvi_crop 5366.1617
julian_day <- gsub("X|_HARV_ndvi_crop", "", row.names(ndvi_df))
julian_day
[1] "005" "037" "085" "133" "181" "197" "213" "229" "245" "261" "277" "293" "309"
origin <- as.Date("2011-01-01")
ndvi_df$date <- origin + as.integer(julian_day)
names(ndvi_df) <- c("tower_ndvi", "date")
ndvi_df
tower_ndvi date
X005_HARV_ndvi_crop 3253.7035 2011-01-06
X037_HARV_ndvi_crop 1888.6754 2011-02-07
X085_HARV_ndvi_crop 2631.4945 2011-03-27
X133_HARV_ndvi_crop 5887.1022 2011-05-14
X181_HARV_ndvi_crop 8793.0962 2011-07-01
X197_HARV_ndvi_crop 8942.3530 2011-07-17
X213_HARV_ndvi_crop 8867.1619 2011-08-02
X229_HARV_ndvi_crop 8856.3461 2011-08-18
X245_HARV_ndvi_crop 8583.9509 2011-09-03
X261_HARV_ndvi_crop 7804.4574 2011-09-19
X277_HARV_ndvi_crop 331.5710 2011-10-05
X293_HARV_ndvi_crop 563.1186 2011-10-21
X309_HARV_ndvi_crop 5366.1617 2011-11-06
ggplot() +
geom_point(data = ndvi_df, aes(x = date, y = tower_ndvi))
field_locations <- read.csv("./data/HARV_PlotLocations.csv")
class(field_locations)
[1] "data.frame"
head(field_locations)
easting northing geodeticDa utmZone plotID stateProvi county domainName
1 731405.3 4713456 WGS84 18N HARV_015 MA Worcester Northeast
2 731934.3 4713415 WGS84 18N HARV_033 MA Worcester Northeast
3 731754.3 4713115 WGS84 18N HARV_034 MA Worcester Northeast
4 731724.3 4713595 WGS84 18N HARV_035 MA Worcester Northeast
5 732125.3 4713846 WGS84 18N HARV_036 MA Worcester Northeast
6 731634.3 4713295 WGS84 18N HARV_037 MA Worcester Northeast
domainID siteID plotType subtype plotSize elevation soilTypeOr plotdim_m
1 D01 HARV distributed basePlot 1600 331.64 Inceptisols 40
2 D01 HARV tower basePlot 1600 341.62 Inceptisols 40
3 D01 HARV tower basePlot 1600 347.61 Inceptisols 40
4 D01 HARV tower basePlot 1600 334.34 Histosols 40
5 D01 HARV tower basePlot 1600 352.93 Inceptisols 40
6 D01 HARV tower basePlot 1600 342.43 Histosols 40
utm18crs <- st_crs(tower)
utm18crs
Coordinate Reference System:
EPSG: 32618
proj4string: "+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs"
field_locations_sf <- st_as_sf(field_locations, coords = c("easting", "northing"),
crs=utm18crs)
field_locations_sf
Simple feature collection with 21 features and 14 fields
geometry type: POINT
dimension: XY
bbox: xmin: 731405.3 ymin: 4712845 xmax: 732275.3 ymax: 4713846
epsg (SRID): 32618
proj4string: +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
First 10 features:
geodeticDa utmZone plotID stateProvi county domainName domainID siteID
1 WGS84 18N HARV_015 MA Worcester Northeast D01 HARV
2 WGS84 18N HARV_033 MA Worcester Northeast D01 HARV
3 WGS84 18N HARV_034 MA Worcester Northeast D01 HARV
4 WGS84 18N HARV_035 MA Worcester Northeast D01 HARV
5 WGS84 18N HARV_036 MA Worcester Northeast D01 HARV
6 WGS84 18N HARV_037 MA Worcester Northeast D01 HARV
7 WGS84 18N HARV_038 MA Worcester Northeast D01 HARV
8 WGS84 18N HARV_039 MA Worcester Northeast D01 HARV
9 WGS84 18N HARV_040 MA Worcester Northeast D01 HARV
10 WGS84 18N HARV_041 MA Worcester Northeast D01 HARV
plotType subtype plotSize elevation soilTypeOr plotdim_m
1 distributed basePlot 1600 331.64 Inceptisols 40
2 tower basePlot 1600 341.62 Inceptisols 40
3 tower basePlot 1600 347.61 Inceptisols 40
4 tower basePlot 1600 334.34 Histosols 40
5 tower basePlot 1600 352.93 Inceptisols 40
6 tower basePlot 1600 342.43 Histosols 40
7 tower basePlot 1600 343.11 Histosols 40
8 tower basePlot 1600 353.55 Inceptisols 40
9 tower basePlot 1600 338.11 Inceptisols 40
10 tower basePlot 1600 343.74 Histosols 40
geometry
1 POINT (731405.3 4713456)
2 POINT (731934.3 4713415)
3 POINT (731754.3 4713115)
4 POINT (731724.3 4713595)
5 POINT (732125.3 4713846)
6 POINT (731634.3 4713295)
7 POINT (731754.3 4713295)
8 POINT (731784.3 4712845)
9 POINT (731904.3 4713595)
10 POINT (731514.3 4713235)
st_write(field_locations_sf, "./data/FieldLocations_HARV.shp",
driver = "ESRI Shapefile", delete_dsn=TRUE)
Deleting source `./data/FieldLocations_HARV.shp' using driver `ESRI Shapefile'
Writing layer `FieldLocations_HARV' to data source `./data/FieldLocations_HARV.shp' using driver `ESRI Shapefile'
Writing 21 features with 14 fields and geometry type Point.